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-^ ABSTRACT 

^^ , Context. It is thought likely that vast numbers of nanoflares are responsible for the corona having a temperature of millions of 

O ,■ degrees. Current observational technologies lack the resolving power to confirm the nanoflare hypothesis. An alternative approach is 

D ' to construct a magnetohydrodynamic coronal loop model that has the ability to predict nanoflare energy distributions. 

^^ ' Aims. This paper presents the initial results generated by a coronal loop model that flares whenever it becomes unstable to an ideal 

MHD kink mode. A feature of the model is that it predicts heating events with a range of sizes, depending on where the instability 
threshold for linear kink modes is encountered. The aims are to calculate the distribution of event energies and to investigate whether 
kink instability can be predicted from a single parameter. 

Methods. The loop is represented as a straight line-tied cylinder The twisting caused by random photospheric motions is captured by 
two parameters, representing the ratio of current density to field strength for specific regions of the loop. Instability onset is mapped 
as a closed boundary in the 2D parameter space. Dissipation of the loop's magnetic energy begins during the nonlinear stage of the 
instability, which develops as a consequence of current sheet reconnection. After flaring, the loop evolves to the state of lowest energy 

ri \ where, in accordance with relaxation theory, the ratio of current to field is constant throughout the loop and helicity is conserved. 

O „ Results. There exists substantial variation in the radial magnetic twist profiles for the loop states along the instability threshold. These 

I . results suggest that instability cannot be predicted by any simple twist-derived property reaching a critical value. The model is applied 

O ' such that the loop undergoes repeated episodes of instability followed by energy-releasing relaxation. Hence, an energy distrilsution 

of the nanoflares produced is collated. This paper also presents the calculated relaxation states and energy releases for all instability 

C/3 , threshold points. 

Cu ■ Conclusions. The final energy distribution features two nanoflare populations that follow different power laws. The power law index 

for the higher energy population is more than sufficient for coronal heating. 

Cn ' Key words. Instabilities — Magnetic fields — Magnetic reconnection — Magnetohydrodynamics (MHD) — Plasmas — Sun; corona 
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fvj ', 1- Introduction ities for how the kinetic energy of these motions is dissipated 

^/-^ - within the loop (Klimchuk 2006). Direct current (DC) heat- 

^ i The Idea that nanoflares are sufficiendy numerous to maintain ■ ^^^^^^ ^^^^ ^^^ ^^^^^^ ^-^^ -^ ^^^^jj compared to the 

in coronal teinperatures was first proposed by Parker (1988). This ^ij^escale of photospheric turbulence. Thus, the loop moves 
O theory implies that the coronal background emission is the result j.^tatically through a series of force-free equilibria until the 

O . of nanoflares occurring continually throughout the solar atmo- accumulating magnetic stresses are dissipated as heat. If the 

• . : 'P"^.^"".^- ".f^^' \^ '' probably unfeasible to observe these flares ^j^^^^^ ^-^^ -^ ^j^^ compared to photospheric motions, alternat- 

individually. This perhaps explains why observational studies ■ ^^^^^^ ^^^^^ ^^^^. ^^^^^ j^^^. fo^t j^t j^otions generate 

have failed to agree on the importance of nanoflares with re- ^^^^^ ^^^^^^^^^ ^^^ ^^^^^^ ^^^^ ^^ ^^-^^ ^ -^ 

L. ■ ^f^l^VT . T?,^ '^f^?!'.? ^^"^ „^?.^^. "^u * ■'''PP within the coronal loop. The waves that can penetrate the corona 

^ , 2000; Aschwanden & Parnell 2002; Parnell 2004). The imprac- ^^ ^^^^^^ ^^ ^^ restricted by type and frequency (Narain & 

■ ■ ' ticahty of nanoflare detection has motivated the development of uimschieder 1996; Hollweg 1984). Whether or not these waves 

inodels that attempt to show how coronal loops might dissipate sufficient energy to heat the corona is not certain (Porter 

their magnetic energy in the form of discrete heating events and ^^ ^j ^994). There is considerably less doubt, however, over the 

thereby replenish coronal heating losses. The primary dissipa- sufficiency of DC heating, which is the basis for the model pre- 

tion mechanism is thought to be magnetic reconnection, since ^^^^^^ ^^^^ ^^^ ^^^^. ^^^jj ^^ ^^^j jble for coro- 

strong evidence of this process has been found in observations ^^^j ^^^^. ^^^ reconnection mechanism mentioned above may 

of large-scale flares (Fletcher 2009; Qiu 2009). Coronal heating ^^^^^^^ ^^^^^ ^-^^^ ^^^ ^^^^^^^ 
requires fast dissipation of magnetic energy, a requirement that 

is compatible with reconnection timescales. The coronal part of the loop is defined by its magnetic 
A loop's excess magnetic energy is introduced via the ran- field (plasma beta, /3 ^ 0.01) and, since the loop's magnetic flux 
dom convective motions that occur at or below the loop's pho- and plasma are frozen together, its magnetic field is continu- 
tospheric boundaries (i.e., footpoints). There are two possibil- ally twisted by random swirling motions at the photosphere. If 
the driving is slow compared to Alfven timescales, the loop's 
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a = {fioj ■ B)/{\B[^) is the ratio of current density to magnetic 
field and r is a position vector (Woltjer 1958). 

Many 3D MHD models have shown how such coronal loops 
exhibit current sheet formation during the nonlinear phase of 
an ideal kink instability (Baty & Heyvaerts 1996; Velli et al. 
1997; Arber et al. 1999; Baty 2000). Essentially, helical current 
sheets become the site of Ohmic dissipation, resulting in a heat- 
ing event. Further simulations have revealed the appropriate cor- 
relation between magnetic energy and Ohmic heating (Browning 
& Van der Linden 2003; Browning et al. 2008; Hood et al. 2009). 
The energy released by an instability depends on the complex 
dynamics of the magnetic reconnection events that occur inside 
the loop. Obviously, this is difficult to model - although this has 
been achieved by 3D MHD simulations. The computational ex- 
pense of these simulations rules them out as a means of explor- 
ing fully the relationship between the o'-profile and the amount 
of energy released. Fortunately, there is another way to calculate 
the energy release. 

Relaxation theory states that when a magnetic field reaches 
instability (or is otherwise disrupted) it will evolve towards a 
minimum energy state such that the total magnetic axial flux and 
the global magnetic helicity are conserved (Taylor 1974, 1986). 
The relaxed state is the well known constant-c or linear force- 
free field: 



VxB = aB. 



(1) 



The original intention of this theory was to explain laboratory 
plasma phenomena; but latterly, it has been frequently applied to 
the solar corona (Heyvaerts & Priest 1984; Browning et al. 1986; 
Vekstein et al. 1993; Zhang & Low 2003; Priest et al. 2005). 
The helicity measures the self-linkage of the magnetic field, see 
Berger (1999). A modified expression for this quantity needs to 
be used since the field lines cross the photospheric boundaries, 
which creates non-zero normal flux and therefore removes gauge 
invariance. 



Jv 



(A +A')(B-B ') dV, 



(2) 



where A is the magnetic potential, B ' is the potential field with 
the same boundary conditions and A' is the corresponding vec- 
tor potential. This relative helicity (Berger & Field 1984; Finn 
& Antonsen 1985) is the difference between the helicity of the 
actual field and that of a potential field (represented by the dash 
terms) with the same normal flux distribution. In ideal MHD, 
the helicity of every flux region is conserved. Taylor proposed 
that in the presence of reconnection all such local invariants 
are destroyed, whilst the global helicity is conserved. Helicity 
is still affected by global resistive diffusion, however, it can still 
be treated as a conserved quantity, if the change in helicity is 
substantially smaller than the change in magnetic energy. The 
results of one aforementioned MHD simulation (Browning et al. 
2008) have shown that 6K/K ~ 10""^ and 5W/W ~ 10"^ (it should 
be noted that these figures are limited by the coarseness of the 
numerical grid used in the simulation). During relaxation, helic- 
ity is redistributed more evenly within the loop; as a consequence 
a becomes invariant with radius. These properties enable one to 
locate a loop's relaxed state and hence calculate the energy re- 
leased during relaxation, i.e., the upper limit of the flare energy. 
In applying relaxation theory to coronal heating, it became 
clear that the effectiveness of the heating depends on how much 
free energy is stored before a relaxation event occurs (Heyvaerts 
& Priest 1984). Browning & Van der Linden (2003) proposed 
that relaxation is triggered by the onset of ideal MHD instabil- 
ity. Thus, the coronal field evolves quasi-statically in response 



to slow photospheric driving, building up free magnetic energy 
- until the field becomes unstable (i.e., the threshold for linear 
instability is reached). At this point, a dynamic heating event 
ensues. It should be noted that ideal instabilities are relevant (as 
opposed to resistive instabilities) because the time-scales are suf- 
ficiently fast. Dissipation and energy release can occur during 
the nonlinear phase of the instability, as has been demonstrated 
extensively by numerical simulations of the nonlinear kink insta- 
bility (Galsgaard & Nordlund 1997; Velli et al. 1997; Lionello et 
al. 1998; Baty 2000; Gerrard et al. 2001). The helical deforma- 
tion of the kink instability generates current sheets in the non- 
linear regime, in which fast magnetic reconnection rapidly dissi- 
pates magnetic energy. Recently, it has been demonstrated using 
3D MHD simulations that this process causes the field to relax 
towards a state which is closely-approximated as a constant-ff 
state, and the energy release is in good agreement with relaxation 
theory (Browning et al. 2008; Hood et al. 2009). These papers 
have investigated field profiles, based on the model of Browning 
& Van der Linden (2003), in the unstable region of parameter 
space and have shown that fast reconnection develops in a cur- 
rent sheet, with dissipation of magnetic energy and relaxation to 
a new, lower-energy state. Efficient "mixing" of the a profile is 
facilitated in the later phases of the evolution, during which the 
current sheet fragments (Hood et al. 2009). 

Supported by this evidence, the model developed here is 
based on the idea that coronal magnetic fields respond to pho- 
tospheric driving by evolving through force-free equilibria, with 
a heating event triggered whenever the ideal instability thresh- 
old is reached. Following Browning & Van der Linden (2003), 
we use a simple cylindrical field model, and represent the nonlin- 
ear force-free field using a two-parameter family of current pro- 
files in which a is a piecewise constant. Browning and Van der 
Linden considered only individual heating events. The primary 
aim of this paper is to predict a distribution of heating events, 
generated by an ongoing stress-relax cycle. The photospheric 
driving stresses the coronal magnetic field until it becomes un- 
stable and relaxes; then the driving resumes, and the process re- 
peats, leading to a series of heating events as expected in the 
nanoflare coronal heating scenario. This is modelled through a 
monte-carlo approach, with the photospheric footpoint motions 
treated as random. Recent observational evidence (Abramenko 
et al. 2006) lends support to the idea that random, turbulent pho- 
tospheric motions can provide an energy source for coronal heat- 
ing in Active Regions. 

In calculating the relaxation and distribution of heating 
events, we need to know the linear instability threshold for line- 
tied coronal loops over a family of field profiles. As a conse- 
quence, we obtain some interesting new results concerning sta- 
bility properties of cylindrical loops, for a much more extensive 
range of current profiles than previously studied (including pro- 
files for which the sense of the twist of the field varies across the 
loop). A secondary aim is thus to explore the properties of the 
linear kink instability on this family of fields, and in particular, 
to determine the extent to which stability can be defined by any 
single quantity such as "critical twist". 

The paper is structured in the following manner. Section 2 
describes the composition of the loop model, along with the 
equations used to express the loop's magnetic field. The calcula- 
tion of the loop's instability threshold is also explained, as well 
as the procedure for allowing the loop to repeatedly undergo in- 
stability and the equations used to determine the energy release 
associated with each relaxation. Section 3 outlines the analysis 
of possible critical parameter values for the onset of loop insta- 
bility. The results of the simulations of energy release are pre- 
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sented in Sect. 4 as nanoflare energy distributions. Finally, in 
the last section, the results are discussed and our conclusions are 
given. 



2. Model 

The model discussed here was used by Browning & Van der 
Linden (2003) and then extended by Browning et al. (2008) to 
include a potential envelope (Fig. [TJ. First, a loop is consid- 




Fig. 1. Schematic of a straightened coronal loop in the r-G plane (left) 
and in the r-z plane (right). The loop, comprises a core (dark grey) 
and an outer layer (light grey); it is embedded in a potential envelope 
(white). The core radius is half the loop radius and 1/6 the envelope 
radius {R1.R2.R3 = 0.5:1:3). The loop's aspect ratio (L/R2) is 20. 

ered to evolve through equilibria as it is driven by photospheric 
footpoint motions. An idealised model of a straight cylindrical 
loop is used with the photosphere represented by two planes 
at z = 0, L; however, the essential physics should apply to more 
complex geometries. The stressed field is line-tied with (in gen- 
eral) a non-uniform a{r) (where a = fiiJ^/B). This is repre- 
sented by a two-parameter family of piecewise-constant-or pro- 
files. Secondly, it is proposed that a relaxation event is triggered 
when the loop's field becomes linearly unstable. The energy re- 
leased, due to fast magnetic reconnection during the nonlinear 
development of the instability, can then be calculated using re- 
laxation theory. 

The MHD kink instability needs to be ideal in order to 
be consistent with the observed rapidity of flare occurrences. 
However, ideal conditions mean there is no resitivity to dissi- 
pate magnetic energy; the 3D MHD simulations discussed in the 
Introduction provide a solution to this impasse. Once a linear 
instability has achieved a positive growth rate, it will soon be- 
come nonlinear: at this point, current sheets will form wherein 
fast reconnection of the magnetic field can take place. These 
expectations are justified by the results of the cited numeri- 
cal simulations. The linear perturbation can be represented as 
/ = fir, zje'^e^', where the azimuthal mode number is set to 1 ; 
this mode has been found to be the least stable (Van der Linden 
& Hood 1999). The effect of such perturbations on the coronal 
loop are represented by the standard set of linearised ideal MHD 
equations. 



2. 1 . Equilibrium fields 

The loop's radial a-profile is approximated by a piecewise- 
constant function featuring two parameters. This design, first 
proposed by Melrose et al. (1994), is readily extensible; extra 
layers of constant a can be inserted to obtain more realistic pro- 
files. The ratio of current to magnetic field is ai in the core, 02 
in the outer layer and zero in the potential envelope. Note that 
the magnetic field is continuous everywhere (though the current 
has discontinuities). Recent work indicates that these a discon- 
tinuities have little discernable effect when compared to similar 
but continuous a-profiles (Hood et al. 2009). 

Without an envelope, the loop's outer surface (located at R2) 
acts as a conducting wall. This is unrealistic in the context of 
the solar corona; the loop would be more stable than it might 
be otherwise. Browning et al. (2008) plotted the relationship be- 
tween the growth rate of the instability and the distance to the 
outer surface of the potential envelope (i.e., a more distant con- 
ducting wall). They found that for six unstable loop states the 
growth rate was invariant once the outer surface of the envelope 
(Rj) exceeded | R2. The R3 boundary was placed at twice this 
value. 

The fields are expressed in terms of the well-known Bessel 
function model, generalised to the concentric layer geome- 
try (Melrose et al. 1994; Browning & Van der Linden 2003; 
Browning et al. 2008). These expressions will change slightly 
whenever o-] or 0-2 become negative; these alterations are cap- 
tured by cr symbols: cri - ^., 0-2 - -r^, and cri2 - cticto. 
(The sign of an a term merely denotes the orientation of the 
azimuthal field). Thus, the field equations for the three regions 
(core, outer layer and potential envelope) are as follows: 



Bi, = BiJo{\ai\r), 

Bw = criBi7i(|ai|r), 

B2: = B2Jo(\a2\r) + C2Yo(\a2\r), 

B20 = (T2(B2Ji(\a2\r) + C2Yi(\a2\r)), 



B'i: — ^3, 

^36 - 0-2 ^2, 

r 



(3) 
Qi<r<Ru (4) 

(5) 
R\<r< R2, (6) 

(7) 
R2<r< R3. (8) 



The fields must be continuous at the inner radial boundaries, Ri 
and R2. Therefore, the constants B2, B3, C2 and C3 can be ex- 
pressed like so: 

^ cTiaMai\Ri)Yo(\a2\Ri) - M\ai\Ri)Yi(\a2\Ri) ,„, 
02 = 01 : , (y) 



^ jQ{\ai\Ri)Ji(\a2\Ri)-cri,2Ji(\ai\Ri)Jo(\a2\Ri) ,,^, 

1-2 = Dl 7 , (lU) 



B3 = B2FQ{\a2\R2), 

C3 = B2^l(|a2l^2), 

where 

2 



7T\a2\Ri 



C2 



Fo.iix) = /o,iW + IT ^0,1 W- 
02 



(11) 
(12) 

(13) 
(14) 
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At all times, the magnetic flux through the loop and envelope is 
conserved: 



r-R?, 

if/* = I Inr* 
Jo 



. . 27rB2 
BUr* ^—^R2Fi(\a2\R2) 
10-2 1 



+ 2KRiBiJi(\ai\Ri)l-^-^ 
\\ai\ ml 

+ nB2Fo{\a2\R2){Rl-Rl\ (15) 

where the asterisks denote dimensionless quantities. Hence, in 
the model, ip* is normalised to 1 and B\ can be determined (not- 
ing that, in equation [TSl B2 is a function of B\). We also nor- 
malise with respect to the coronal loop radius, R2, see Fig.[T] 

As the random motions of the photosphere proceed, the loop 
evolves through a series of force-free equilibrium states until it 
becomes linearly unstable. We now discuss the calculation of the 
instability onset. 



2.2. Linear kink instability ttireshold 

A coronal loop's instability is constrained by the line-tying of 
the photospheric footpoints (Hood 1992). Hence, all perturba- 
tions are required to vanish at the loop ends (z = Q,L). When 
the growth rate of a perturbation transitions from a negative 
value to a positive one, the loop has reached the threshold of 
an ideal linear instability. The instability threshold is a curve 
in 2-dimensional a-space (ai, 02). The properties of the loop 
(e.g., ai and 02) at these threshold points can be found by sub- 
stituting the perturbation function into the linearised MHD equa- 
tions, leading to an eigenvalue equation for the growth rates 
(Chap. 7, Priest 1987). The growth rates and eigenfunctions of 
the most unstable modes are found numerically, for line-tied 
fields, with the CILTS code, described in Browning & Van der 
Linden (2003) and Browning et al. (2008). CILTS can be config- 
ured such that one of the loop's a parameters is fixed whilst the 
other is incremented. The code terminates as soon as the real part 
of the eigenfunction falls below zero, i.e., the loop is no longer 
unstable to kink perturbations. The left panel of Fig.|2]shows the 
closed instability threshold curve mapped by the CILTS code 
(see also Fig. 5 of Browning et al. 2008). The threshold curve 
has symmetry: it is invariant when rotated by n radians. Thus, 
it is sufficient to show how various properties (e.g., magnetic 
twist and energy release) vary along the top half of the thresh- 
old curve. For ease of plotting we can convert this half of the 
threshold curve to a one dimensional form: the filled circles and 
bold numbers shown in the right panel of Fig. |2]represent the tic 
marks and labels for the ID threshold point axis, see Figs.lSl-fTTl 
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Fig. 2. The left panel shows the closed instability threshold (solid) with 
the B. reversal lines (dashed). The top half of the threshold (where 02 > 
0), annotated with threshold point numbers, is shown in the right panel. 

It is important to remember that the threshold only applies 
to the specific loop geometry outlined above. A new threshold 
would need to be calculated should the loop's proportions, com- 
position or envelope change as a consequence of some activ- 



ity. Another caveat is that there are some points in a-space that 
yield singularities when calculating quantities such as helicity 
or magnetic energy. These arise when the axial field (B.) has a 
significant region of reversal; since i//* is normalised to 1 and 
conserved, the helicity and energy terms will diverge when un- 
normalised i/r* — > (see Browning & Van der Linden 2003). 
Fortunately, the instability threshold does not enter the region 
where B. begins to pass through zero (see dashed lines in Fig. 
|2]i, so these singularities are not encountered. 

2.3. Random walk 

When a loop is twisted by turbulent photospheric motions, it 
performs a random walk through the a-space enclosed by the 
instability threshold (Fig. [3]). This traversal of cc-space is ran- 
dom in direction but constant in length (we set the dimension- 
less step-length to be 6a = 0.1 for most of the results presented 
here). Clearly, the nature of this random walk depends on the 
statistical properties of the driving photospheric motions; in fu- 
ture, different forms of this random driving will be investigated, 
but for now, we take the simplest assumptions. 

The time unit r is the step time, the time taken for a to change 
by Q.I/R2 (in dimensional units). We may estimate, roughly, 
a timescale for this process as follows. Based on axial val- 
ues, a change 6a corresponds to a change in magnetic twist 
d(p ~ {L/2)(6a/R2y, taking L/R2 = 20 gives 50 ~ 1. If this is caused 
by photospheric twisting motions of magnitude vg for a time in- 
terval T, we find T ^ {6(p)Rf/vg, where Rf is the footpoint radius. 
With typical values of Rj = 200 km and ve = 1 km s"' , we obtain 
T ^ 200 s; note that this is consistent with quasi-static evolution, 
justifying a posteriori our choice of random-walk step size. The 
step time (r) may be identified with the correlation time of pho- 
tospheric motions, which is likely to be rather longer than the 
value given above; for example, a granule lifetime of 1000 s may 
be appropriate (Zirker & Cleveland 1993). The effect of increas- 
ing the step length (da) is considered in Sect. 14.2.31 

Eventually, the field will reach the instability threshold: it 
will become linearly unstable. At this point, the field releases 
energy and transitions to a lower-energy state defined by Taylor 
relaxation: helicity is conserved and the cf-profile relaxes to a 
single value. 



2.4. Energy release calculation 

We have extended the model (Browning & Van der Linden 2003; 
Browning et al. 2008) to allow a loop to repeatedly undergo re- 
laxation as it evolves within a-space. Initially, a loop starts from 
a randomly-selected stable state. The field profile then undergoes 
a random walk until it crosses the instability threshold; where- 
upon, the loop relaxes and the profile transitions to the relax- 
ation line (a\ = 02). The constant a-value (q'„) will, of course, 
vary depending on where the threshold was crossed; q;„ is found 
by helicity conservation (Browning & Van der Linden 2003). In 
mathematical terms, we find the roots of the following equation: 

K{a,.,)-K{ai„,ai,2) = 0, (16) 

where am and a,72 are the coordinates of the instability thresh- 
old crossing (conservation of axial flux is assured through the 
normalisation (/»* = 1). The helicity can be expressed as follows: 



K = 2L 



nRi 
Jo 



mr) 



dr. 



(17) 



where / is the current and L is the loop length (Finn & Antonsen 
1985). The equation for the magnetic energy contained within 
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the loop and envelope is straightfoward; 



Ln r*' , 
W = — rB-dr, 

y"0 Jo 



(18) 



where L is normalised to 20 (since Ro — V) and //q is set to 1. 
See Appendix A for the full expressions. The energy difference 
between the unstable and relaxed states can be calculated thus: 



6W = W(ai,i,aia)-W{a,.,). 



(19) 



This is the relaxation energy: the energy released as heat during 
the event. 

After relaxation, the loop resumes its random walk until it 
reaches the threshold and the process repeats. Fig.[3]gives an il- 
lustration of this process. We thus generate a sequence of energy 




Fig. 3. The instabihty threshold encloses the relaxation line, which is 
a subsection of the a\ = ai line centred on the origin. The annotations 
illustrate the initial stages of a simulation that begins at the position 
marked by the cross. The first random walk is shown in light grey; it 
ends at the threshold position marked ▲ and the associated relaxation 
point is indicated by A, which is the starting point for a second ran- 
dom walk (dark grey). This walk attains instability onset at the position 
marked ■ and relaxes to the point labelled n - the starting point for a 
third walk. 

release events, which are collated to produce a nanoflare energy 
distribution, see Sect. 14.21 

When the loop relaxes, the a-profile throughout the loop and 
envelope becomes constant, i.e., a\ —a2-a-ii^Q- The envelope 
is no longer potential; it has acquired a residual current. In prin- 
ciple, as the main portion of the loop is twisted again by on- 
going motions, a new equilibrium would develop with varying 
current {a\,a2) in the loop and a non-zero current {aj) in the 
envelope. For some threshold sections the consequent residual 
current is so small that the threshold shape would remain un- 
changed. However, the validity of the simulation process can 
only be assured by including an extra stage, wherein the enve- 
lope dissipates its helicity so that it becomes potential again. The 
loop can now resume its random walk with respect to the same 
threshold. 

Although the primary purpose of this model is to calculate 
the distribution of energy releases, in achieving this we also ob- 
tain some interesting new results on linear stability, which are 
summarised in the next section. This is because, in order to ex- 
plore the full parameter space of equilibrium current profiles, we 
have calculated the linear stability properties of a much wider 
family of fields than previously investigated: in particular, fields 
with reversed twists. 



3. Instability threshold and critical twist 

The evolution of the field profile through a-space is determined 
by photospheric perturbations which map to changes in the mag- 
netic field and hence to changes in a\ and q'2- A loop's magnetic 
twist is directly related to rotational photospheric motions: 



V 



rB,{r) ' 



0<r<Ri. 



(20) 



Thus, the photospheric motions directly determine a loop's ip- 
profile, which in turn determines a{r) (and hence, in our model, 
ffi and Q'2). The magnetic twist of coronal loops is an observable 
feature (Kwon & Chae 2008). Portier-Fozzani et al. (2001) have 
even observed a loop's twist decreasing over time - evidence per- 
haps of a loop evolving towards a state of minimum energy. 

3.1. Criteria for instabWity 

Many workers have looked at the idea that the magnetic twist of 
a loop can be used as a proxy for the onset of kink instability with 
some critical twist parameter determining instability. Hood & 
Priest (1979) performed a MHD stability analysis on a variety of 
straightened line-tied coronal loops. They showed that the criti- 
cal twist, (fcrii, varies according to the aspect ratio {L/R2), plasma 
beta (J3) and the transverse magnetic structure (i.e., the nature of 
the variable twist profile). In general, short fat loops (low aspect 
ratio) have low critical twists, whereas long thin loops (high as- 
pect ratio) have high critical twists. Subsequent work (Hood & 
Priest 1981) revealed that a loop of uniform twist, aspect ratio 1 
and zero j3 possessed a critical twist of 2.497r. This figure is often 
quoted as a general result. 

The question arises as to whether there is any single param- 
eter (such as peak or average twist) which determines instability 
onset for all twist profiles. Indeed, more generally, it would be 
desirable to have a single quantity determining instability onset 
even for more complex (non-cylinderical) fields (Malanushenko 
et al. 2009). Incidental to our task, we have calculated the sta- 
bility properties of an extensive family of equilibria, including 
fields with reversed twist (as well as simple monotonic-twist pro- 
files as used by other authors). This provides a very useful test 
for any proposed criteria for instability onset. 

Loops with variable-twist profiles have been previously stud- 
ied; however, all such profiles generally have a similar form 
(Vein et al. 1990; Mikic et al. 1990; Baty 2001): the axial twist 
is the maximum, then the twist declines to a negligible value 
at the loop boundary. Velli et al. (1990) calculated that insta- 
bility occurred when ipo = (p{r=Q) = 2.57r; this agrees with Hood 
and Priest's result for a uniform twist profile. Mikic et al. (1990) 
calculated a critical axial twist of 4.87r, the loop's average twist 
however, was ~2.57r. 

The idea of using magnetic twist as a marker for instabil- 
ity relies on the existence of some twist-derived parameter hav- 
ing a constant value for all the points on the instability thresh- 
old. Baty (2001) used a MHD stabiHty code to show that for 
a small set of equilibria the average twist at instability is the 
same for several different magnetic configurations. However, 
there are several differences between Baty's work and the model 
presented here. Firstly, the twist profiles defined for each equi- 
librium are all positive and none contain multiple peaks (Fig.l, 
Baty 2001). Furthermore, critical twist convergence arises when 
the normalised distance, d, is greater than 5 (Fig. 5, Baty 2001). 



d = 



V0^2 



(21) 
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In our equilibria, (i = (/3o/20, the threshold values of the absolute 
axial magnetic twist vary from zero to 9n, which means d varies 
from to 1 .42, thus the average twists will not be in the regime 
where the critical axial twist converges to I.Stt. In fact, the large 
d regime cannot be attained. 

The idea of a single critical average twist seems unlikely 
when one examines the threshold presented here. Clearly, at 
some threshold points a^ and 02 are of opposite sign; thus, one 
or two places on the threshold will have an average twist equal to 
zero, and yet they are unstable. Perhaps critical twist is achieved 
within a subsection of the loop; this idea is explored further in 
Sects.lOandllS] 



3.2. Radial twist profiles and linear eigenfunctions 

In order to understand the nature of the instability, we investi- 
gate the twist profiles and eigenfunctions of the unstable mode, 
for different parts of the instability threshold. The magnetic twist 
profiles for a selection of points just outside the threshold curve 
exhibit considerable variation, as Fig. |4] illustrates. For the un- 
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Fig. 4. The loop's radial twist profile at specific points along the thresh- 
old (labelled A-F). 

stable equilibrium labelled as point B (ffi < 0, Qf2 > 0), the cor- 
responding twist profile (also labelled B) shows that the core 
field has a strong negative twist, whilst in most of the outer layer 
and in all of the envelope the field has positive twist. As one 
moves from A to B, the twist in the core becomes more neg- 
ative whilst the twist in the outer layer moves in the opposite 
direction. It appears that the increase in a2 stabilises the neg- 
ative core twist by providing additional reversed (i.e., positive) 
twist in the outer layer The sharp corner at the top left of the 
threshold marks the point where instabilities driven within the 
core intersect those that originate from within the outer layer 
Profiles C to D therefore, suggest instabilities driven in the outer 
layer, since \a2\ > \ai |. The peak twist in the outer layer reduces 
as the core twist moves from negative to positive. Further along 
the threshold, where ai > q'2, the instabilities are likely to be 
driven in the core. Note that profiles D, E and F are always posi- 
tive in sign; D has a twist peak near the loop edge while E and F 
are roughly monotonically decreasing. The corresponding mag- 
netic field profiles for the six points A-F are given in Appendix 
B. 

It seems that instabilities are driven mainly on or near the 
peak of largest absolute twist. A twisted field region may be sta- 
bilised by an enclosing region of opposite twist. Furthermore, 
a twisted outer layer may need less twist to achieve instability 



if the core has the same twist orientation. To investigate these 
ideas further, we plot the unstable eigenfunctions obtained from 
CILTS for the same a-space points, A to F. 





Fig. 5. The linear eigenfunction, Vt(x,y=0,z), for the a-space points 
A (left) and B (right) profiled in Fig.|4] Cartesian coordinates are used, 
hence, the x-axis is equivalent to the radial axis. 

The eigenfunctions for profiles A and B (Fig. |5]l show that 
the amplitude is strongest in the core. Interestingly, the ampli- 
tude for profile A has dropped to zero long before the envelope 
boundary at R2, which suggests that in the subsequent relaxation 
only inner regions will be affected, with little change in the po- 
tential envelope. There is a strong similarity between the eigen- 





Fig. 6. The linear eigenfunction, V^(x,y=0, z), for the ff-space points C 
(left) and D (right) profiled in Fig.H] 

functions for profiles C to D (Fig.|6]l and the amplitude is high- 
est near the R2 boundary, indicating an outer layer instability. 
Finally, the two plots in Fig. [T] (profiles E and F), clearly show 
a progression towards the eigenfunction calculated for A (albeit 
with a Vv of opposite sign). Notice also that the form of the 





Fig. 7. The linear eigenfunction, V^.(x,y=0, z), for the a-space points E 
(left) and F (right) profiled in Fig.|4] 

eigenfunction changes significantly between B and C, indicating 
that diff'erent modes are going unstable. This is to be expected, 
since the ff-space positions labelled B and C (Fig. H] top left) 
are either side of the intersection point formed by the two curves 
that describe the instability threshold. 
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3.3. Critical twist parameters 

Next, we look for a twist-related parameter that takes on a criti- 
cal value whenever the loop reaches the threshold. As expected, 
the variation in axial twist, i^o, is similar to the variation in ip{Ri), 
see Fig. [8] None of the quantities suggest any single (constant) 
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Fig. 8. The variation in magnetic twist around the instability threshold 
(the threshold points are defined in Fig. [2} for three radial positions. 
The solid line represents the variation in axial twist, ipg; the dashed line 
is the variation in twist at the boundary between the core and the outer 
layer, tp(R[); and the long-short dashed line is the variation in twist at 
the boundary between the outer layer and the potential envelope, (fiiRi). 
The twist values are plotted in units of n. 

critical value. Perhaps, the average twist is less variable around 
the threshold? There are several ways to calculate this expres- 
sion; 



{^)q 






<'^>o' 



~ nR^ Jo 



rBSr) 



2nr -— dr 

rB~Sr) 



(22) 

(23) 
(24) 



where Ri is a radial upper bound (e.g., R2 or R^). Equationl24lis 
the average twist weighted by area. The other two equations (|22| 
and|23]) have been used by Baty (2001) and Velli et al. (1990). 
Note, Equation |22] can be calculated analytically, see Appendix 
A. {tp)^ denotes the average twist, weighted by area, over the 

core and outer layer Similarly, (^)o' denotes the same quantity 
but over the loop and potential envelope. The tilde (~) and hat 
(A) symbols are used to indicate the other equations. 

When ffi and ai are equal there is no distinction between the 
loop regions; this occurs on the threshold when ai = a2 ~ 1 -4 - at 
this point {ip)^ ~ 5k. When 0-2 = and ai * 2.2 the loop is identi- 
cal to the core with a bigger potential envelope and (^)q' ~ 7.7 n. 
Thus, as expected, fatter loops like the first case, have lower in- 
stability twists than thinner ones. When the loop's core is poten- 
tial (i.e., when ai-0 and a2 ~2.2), (v')^' a;4.57r is the average 
outer layer twist. In this configuration, the loop is less stable than 
the case when only the core is non-potential. 

None of the twist averages (Fig. |9]l is invariant along the en- 
tire threshold. Although, {ip}q^ and (ijo),,' have approximately the 
same value {^ l.ln) between threshold points 40 and 90. All 
of these points show a positive peak twist at R2. Hence, the en- 
velope's contribution dominates the overall average twist. This 
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Fig. 9. The variation in the average twist over the core, 0-R[ (top left), 
the outer layer, R\ -R2 (top right), the loop, O-R2 (bottom left) and the 
loop and envelope, - R3 (bottom right). The solid lines were calculated 
according to Eqn. |24l the dashed according to Eqn. |23] and the long- 
short dashed according to Eqn. [22] Again, the twist values are plotted 
in units of n. 



is especially true for the (v3)q' case: the higher the radial coor- 
dinate the greater the weight of the calculated twist. Within the 
envelope, the twist declines as 1/r and so, the inclusion of the 
envelope twist averages out the final result. 

Notice also, that all the twist averages go through zero when 
the core and outer layer have opposite twists. It seems that these 
quantities do not reveal the detail necessary to understand why 
a particular loop configuration is on the point of instability, nor 
where in the loop that instability originates. 

Finally, we consider the proposal of Malanushenko et al. 
(2009), that a critical value of normalised helicity (equivalent, in 
our terms, to the normalised loop helicity, K/ifP', over the range 
- R2) indicates instability onset. In fact, the normalised helicity 
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Fig. 10. The variation in Kji^f^ (over the range O-R2) along the insta- 
bility threshold. The threshold states between the vertical dashed lines 
feature reverse twist; outside the lines the twist is single-signed. 

is certainly not the same for every threshold point; this quantity 
passes through zero because ai and 02 take on values of oppo- 
site sign along some sections of the threshold. For fields with 
single-signed twist, the normalised helicity gives an approxi- 
mate threshold, but even here, the (absolute) critical value ranges 
from about 1 .5 - 2.5. Figure [TO] shows that the idea of such a crit- 
ical value breaks down for loops that feature regions of reversed 
twist. 
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4. Distribution of energies and coronal heating 
considerations 

We now proceed to the main task of our work, which is to cal- 
culate the distribution of magnitudes of the sequence of heating 
events generated by random photospheric driving. 

4. 1 . Helicity and energy 

The top left panel of Fig. [TT| plots total helicities of the thresh- 
old states. A total helicity (or flux) is one calculated over the 
range O-R},, i.e., the loop and envelope. None of the threshold 
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Fig. 11. Helicity (top left), magnetic energy (top right), relaxed alpha 
(bottom left) and energy release (bottom right) along the ID representa- 
tion of the instability threshold. The energies, W* and 6W', are dimen- 
sionless quantities. 

states have sufficient helicity for the relaxed state to feature he- 
lical modes (Taylor 1986) and so all relaxed states are cylin- 
derically symmetric. The bottom left panel confirms that each 
threshold state corresponds to a relaxed state. The maximum »„ 
is about 0.48. Hence, there is a good chance that the envelope 
current after relaxation will not be insignificant. This result con- 
firms the need to investigate how this current could be dissipated 
(although the maximum value is still significantly less than typ- 
ical ffi and Q'2 threshold values). 

The energies shown in Fig. [TT] (top right and bottom right) 
are given as dimensionless quantities, to calculate the dimen- 
sional energy we first give the dimensional flux though the loop 
and envelope. 



^ = 



3R,. 



2nrB-dr = 9nRlBr 



(25) 



where Re is the coronal loop radius (the dimensionahsed Ro) and 
Be is the mean axial coronal field. Since we have normalised the 
magnetic field, such that kjj* -I, the field strength can be dimen- 
sionalised thus. 



B = -^B' = 97iBcB\ 
Rl 



(26) 



where the asterisk is again used to denote dimensionless values. 
Hence the dimensional energy release becomes 



6W 



111. 



R' dW* 






rIbI 6W* . 



(27) 



This expression differs slightly from the one used by Browning 
& Van der Linden (2003): their dimensionless energy release is 
calculated per unit length and Rj, -R2- Assuming typical values 
{Re - 1 Mm and Be = 0.01 T), we obtain dimensional energy val- 
ues of 6 X 10^2 6W J = 6 X 10-'^ 6W erg. Thus, the top end of 
the SW* scale («; 0.073) is equivalent to 4 x 10^^ erg. This is in 
the microflare range, but nanoflare energies will be obtained for 
weaker fields or for smaller loops. 

4.2. Flare energy distribution 

Every time a loop's random walk reaches the instability thresh- 
old, the loop's relaxed state (i.e., its position on the relaxation 
line) and energy release are calculated. The relaxed state is the 
start of a new random walk which will lead to another relaxation 
(and energy release). If this process is repeated often enough it 
can be shown that certain energy release sizes are more common 
than others. The flare energy distributions converge as the num- 
ber of relaxation events simulated (see Sect. 12.4b is increased. 
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Fig. 12. Flare energy distributions for 100 (top left), 1000 (top right), 
lO** (bottom left) and 10^ (bottom right) relaxation events. 

The gross features of the converged energy distribution can 
be explained by presenting the energy distribution for the high- 
est number of flare events (Fig. [12] bottom right) alongside the 
instability threshold in Fig. [13] Both the distribution and the 
threshold are colour-coded according to event energy. The colour 
of the threshold point encountered by a coronal loop indicates 
the energy of the resulting flare and also the part of the distri- 
bution where the heating event falls. Dimensionalised values for 
the flare energies have been recovered by using typical loop pa- 
rameters (see above). 

The energy release changes as one moves along the thresh- 
old, i.e., there is an energy release gradient. This gradient is 
small for low energies, therefore, only a few bins cover the low 
energy sections of the threshold, see Fig. [13] The low energy re- 
lease events are divided amongst a small number of bins, hence, 
these bins contain many more events. The threshold sections cor- 
responding to the profile minima are slightly longer than those 
associated with the first peak [^ 1.5 times), however, the minima 
sections have a much higher energy release gradient. The energy 
release events are divided amongst a higher number of bins and 
so each of these bins contain fewer events. As one moves into 
the threshold sections that correspond with the second profile 
peak (denoted by blue-green shades) the energy release gradient 
decreases. Hence, one would expect the associated bins to have 
higher event counts. This increase is accentuated however by the 
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proximity of the corresponding relaxation points. Once a loop 
achieves a blue-green instability, it will relax to a point (also 
coloured blue-green), which happens to be closest to the green 
section of the threshold. Subsequent walks will have a higher 
chance of crossing that section, thus, the corresponding distri- 
bution bins will contain even more events. The highest energy 
releases available on the threshold are farthest away from the 
relaxation line. In this part of the distribution, the chance of an 
energy release event is inversely proportional to the size of the 
energy release. 
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Fig. 13. The flare energy distribution for 10' relaxation events and the 
instabihty threshold (with relaxation line). All lines are colour-coded 
according to energy release. 
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Fig. 14. Flare energy distributions over 10' relaxation events for a vari- 
ety of loop lifetimes. The lifetimes are 1000 relaxation events (top left), 
100 relaxation events (top right), 10 relaxation events (bottom left) and 
1 relaxation event (bottom right). 

These results are not strongly tied to the loop lifetime, which 
is the number of relaxations a loop undergoes during the simu- 
lation. We can simulate loop replacement by randomly selecting 
a new position within the stability region after a certain number 



of relaxations. The top left distribution of Fig. [14] is generated 
from a loop that is replaced after every 1000 relaxation events. 
As the loop lifetime is reduced so is the height of the second 
peak. However, this peak reduction is only noticeable for small 
lifetimes, e.g., < 10 relaxations. A minimum lifetime of 1 relax- 
ation (Fig. [141 bottom right) still yields a two-peaked distribu- 
tion, albeit with a reduced second peak. This result corresponds 
to finding the energy release from an ensemble of identical loops. 
It is also the most conservative in terms of total energy released. 
The concept of a loop ensemble (a collection of 10^ loops 
flaring simultaneously) is particularly useful, since it allows us 
to sidestep the complications that come with allowing loops to 
survive many relaxations. Otherwise, we would need to under- 
stand how a loop is affected by its energy release. A loop may 
shrink or implode after flaring (Janse & Low 2007). The instabil- 
ity threshold would be invalidated should the loop's aspect ratio 
be altered as a consequence of post-flare shrinkage. Therefore, 
the next section will examine in more detail the ensemble distri- 
bution. 



4.2.1 . Distribution of "nanoflares" 

The ensemble distribution (Fig. [141 bottom right) does not yield 
a simple inverse power-law when converted to a log scale; there 
are two peaks in the profile. The trailing edge of the first peak 
equates to a power-law slope of x; -1.5. Its internal structure can- 
not currently be resolved since one would need to increase the 
threshold point density. The trailing edge of the second peak 
gives a slope of ^ -8.3; this is much greater than the critical gra- 
dient for nanoflare heating, m < -2 (Hudson 1991). 

These power-law figures are provisional and are likely to 
change as the model is enhanced. For example, a more realistic 
or-space traversal function - one where 6a2 is correlated with 6ai 
- will prefer walks parallel to the relaxation line; this will alter 
the distribution of heating events. Furthermore, we consider here 
only a single loop of fixed dimensions (or an ensemble of iden- 
tical loops). In reality, in any given large-scale loop structure, 
sub-loops of varying radii will be generated, depending on the 
horizontal scale-length of the driving photospheric motions. The 
distribution must then be averaged over a distribution of radii, as 
well as considering variations in length and field strength. 

4.2.2. Heating Flux 

The primary aim of our model is to calculate the distribution 
of nanoflares. In general terms, the heating rate will be similar 
to other calculations in the literature based on random photo- 
spheric twisting (Abramenko et al. 2006; Zirker & Cleveland 
1993; Berger 1991; Sturrock & Uchida 1981). Within our ap- 
proach, all the energy input from the photosphere must be dissi- 
pated, in a long-term time-average over many events, since the 
build up of coronal magnetic field is limited by the instability 
threshold. 

The energy flux, F, can be expressed using Eqn.l27l 



F = 



!i![! ± ^ ^3^? ^,^.^ ^ 81. ^ (^ 



where A^ is the average number of steps taken to reach the thresh- 
old, r is the time taken to complete each step in the random walk 
and {6W*) is the average dimensionless energy release, given by 
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For the ensemble case, {6W*}-- 
Applying previously used values (Be ■ 



Rc = l Mm), yields a total flux of 6 x 10 erg cm ^ s 



10.0293 and A^«264. 

:0.01 T, T = 200 s and 
This 

result is applicable to Active Regions (a value for the Quiet 
Sun can be obtained by setting Be = 0.001 T; this simply low- 
ers F to 6 X 10"* erg cm - s '). These results are slightly lower 
than Withbroe & Noyes (1977) measurements of coronal heating 
losses, which are 10^ erg cm"^ s"' (Active Regions) and 3x10^ 
erg cm"^ s ' (Quiet Sun). This shortfall vanishes however if the 
random walk is conducted using larger step sizes, see following 
section. 



4.2.3. Random walk step size 

So far, we have assumed photospheric motions to be somewhat 
temporally incoherent, with a short random walk step, 5a = 0.1, 
corresponding to t = 200 s. Here, we consider the effect of vary- 
ing this step length; in other words, varying the coherence time 
of the photospheric motions. Fig. [15] which should be compared 
with the bottom right panel of Fig. [141 shows the distributions 
that result when Sa is increased by factors of 10 and 40. For 
the latter case, the step is sufficiently long that, on average, just 
one step is required to reach the instability threshold: this corre- 
sponds to a temporally coherent twisting of the loop (although 
the spatial profile of the twisting varies randomly between heat- 
ing events). 
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Fig. 15. Flare energy distributions over 10^ relaxation events for a step 
size of 6a = I (left) and Sa = 4 (right). Loop lifetime is one relaxation 
event. 

It can be seen that the energy distribution is virtually inde- 
pendent of the random walk step size. However, the heating flux 
does depend on this quantity. It is expected that the heating flux 
should be proportional to r (Berger 1991; Zirker & Cleveland 
1993). Eqn.l28]shows this to be the case, since the average num- 
ber of steps, A^, for the random walk to reach the threshold scales 
as A^ oc T"2. In particular, if we consider the limiting case of co- 
herent twisting (t = 8000 s, N^l.U and <5W*)«0.0305) we 
find that the flux increases significantly, F^3xlQ^ erg cm"^ 
s"' for Active Regions. 



4.2.4. Random walks in twist space 

The instability threshold can be expressed in terms of ip{R\) and 
(piRi)- This is more representative of the twist profile (see Sect. 
[3]l, which is directly generated by photospheric motions. We use 
(f{R\) and ip{R2) to obtain a 2-parameter representation of the 
family (p{r), since the twist profiles of Fig.|4]usually show peak 
twists at the radial boundaries/?! (0.5) and/?2 (1-0). The instabil- 
ity threshold can thus be plotted in twist space rather than alpha 
space, as in Fig. [3] Given that the process is driven by turbu- 
lent photospheric motions, it is more realistic to consider that 
the twist randomly evolves, rather than the a profile. We thus 
repeat our calculations with a random walk in i^-space. 
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Fig. 16. The instability threshold and relaxation line in ^-space (left), 
alongside the flare energy distribution for a 10' loop ensemble per- 
formed within i/3-space (right). The random walk step size, 6(p, is ap- 
proximately 0.32;r, this corresponds to a step time of 200 s. 

When a simulation is run, for a sequence of heating events, 
in ^-space, the resulting energy profile is more or less identical 
to that generated by an a-space simulation, see Fig. [T6l (right). 
The translation to ip-space results in a slight reorientation of the 
relaxation line with respect to the threshold. Consequently, the 
relaxation line is closer to the part of the threshold associated 
with the highest energies; this explains why the distribution has 
a thicker tail. The energy flux is F x4x 10^ erg cm"^ s ', which 
increases with step size in the same manner as shown for the 
Qf-space simulations. 



4.3. Temporal properties 

We can gain some insight into the temporal distribution of heat- 
ing events by assuming that each step of a random walk in a- 
space takes the same arbitrary unit of time, corresponding to 
photospheric driving. The time axes of Fig. [17] are shown in 
terms of random walk step number (i.e., step count or running 
step count). The time unit r is the time taken for one step, which 
was estimated to be 200 s, see Sect. 12.31 
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Fig. 17. Top Left: the number of steps between relaxation and insta- 
bility (i.e., time interval between heating events) for each event of a 
simulation comprising lO'* events. Top Right: the number of steps taken 
to reach the threshold against the energy released. Bottom Left: the en- 
ergy release as a function of time (running step count). Bottom Right: 
the 10' event simulation produced 10' flares of varying energies. The 
size of the flares are shown in a way that is reminiscent of actual 
flare/microflare/nanoflare observations: the bigger the event, the wider 
the base of the triangle used to represent that flare. The figure covers a 
time sequence equal to 5000 steps, taken from a random position within 
the simulation data. 

The probability that a flare event will have occurred after 
a particular number of steps in ff-space is invariant with time, 
see Fig. [17] (top left). This is also true for the probability that a 
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flare event wifl have a particular energy: the flare energy is not 
dependent on where in the simulation it occurs. The two hori- 
zontal bands shown in the bottom left panel are consistent with 
the peaks shown in Figs. [T4lfT5] 

There does appear to be a slight relationship between the 
flare interval time (i.e., the number of steps taken to reach the 
threshold) and the flare energy, see Fig. [17] (top right); but, no 
positive correlation is evident. Observations strongly suggest 
that these two properties are uncorrelated (Wheatland 2000). 
The top right panel echoes the shape of the energy release distri- 
butions. The instability threshold has an energy release gradient, 
which means that sections of threshold that have a more or less 
constant energy release (low gradient) will be visited by more 
loops. The higher the number of visiting loops, the wider the 
range of step counts associated with that section of the thresh- 
old. Thus, Fig.[T7](top right) is in agreement with observations. 

Finally, the simulated energy releases can also be represented 
as a time series of flare energies, again see Fig. [TTl (bottom left). 

4.4. Critical magnetic sliear and coronal heating 
requirements 

Parker (1988) explains that a magnetic flux bundle rooted in 
the photosphere (i.e., a coronal loop) will be stressed as pho- 
tospheric motions move the footpoints in a direction transverse 
to the original magnetic field. The coronal heating requirements 
for Active Regions and for the Quiet Sun can be converted to 
equivalent magnetic stresses, which can in turn be converted to 
magnetic shears {B±/B~); 

B^ 4nF 

(30) 



B 



B2 



V0 



where B^ is the transverse field, B^ is the original axial field, F is 
the coronal heating requirement and vg is the photospheric flow 
velocity. We set vg = 1 km s"' for consistency and find that the 
necessary magnetic shear is approximately 0.12 (~7°) for Active 
Regions and 0.38 (-21°) for the Quiet Sun. 

Thus, for coronal heating to work, there must be some mech- 
anism which restricts the shear to around these levels. If dis- 
sipation occurs at lower levels of shear the energy flux cannot 
maintain coronal temperatures. Conversely, if the shear can build 
up to much larger levels, the energy input would be higher than 
required. Our model provides an explanation for these critical 
shear values. We calculate two types of mean magnetic shear 
along the threshold, the mean absolute shear and the root mean 
square shear: 



-AT 



dr , 



\ 



2 r«^ lBe\~ 



dr . 



(31) 



(32) 



Both quantities are area-averaged. The plot (Fig.fTSll shows that 
the values of Smabs and Srms are comparable with those de- 
rived above. In general, we find that the average shear at the 
threshold is slightly higher than the limit suggested by Parker 
Perhaps more shear is required, since not all of the magnetic en- 
ergy released during relaxation will be converted to heat. Fig.fTSJ 
suggests that the apparent limiting value for magnetic shear is 
determined by the linear instability. 

We should mention the work of Dahlburg et al. (2009), who 
provide an alternative explanation of this shear dependency in- 
volving an explosive "secondary instability". They have used a 
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Fig. 18. The mean absolute (soUd) and root mean square (dashed) of 
the magnetic shear along the instability threshold. The shears are calcu- 
lated over the loop volume, Q-R2. 

3D viscoresistive MHD code to simulate the shearing of a line- 
tied flux tube. The results of this code have revealed that the 
secondary instability requires a critical shear for onset. 

5. Summary and conclusions 

We have used a model based on relaxation theory to determine 
the energy release of a sequence of heating events generated 
by random photospheric footpoint motions. This is a "stress- 
relax" scenario, in which free energy is continually built-up in 
the corona due to photospheric driving - repeated relaxations oc- 
cur whenever a critical condition is reached. We postulate that 
heating events are triggered by the onset of an ideal kink insta- 
bility. Energy release then occurs by fast magnetic reconnection 
in current sheets generated during the nonlinear phase of the in- 
stability, and this is calculated by assuming a helicity-conserving 
relaxation to a minimum energy state. The energy release of any 
individual heating event depends on the current profile at the 
point where the instability threshold is crossed - thus heating 
events of a range of sizes are found. 

The initial stressed fields are modelled using piecewise- 
constant profiles of a, yielding a two-parameter family of cur- 
rent profiles. Photospheric driving thus yields a random walk of 
the current profile. In order to pursue the calculation, the ideal 
instability threshold of a line-tied loop has been determined. 
The modelled current profiles incorporate regions of positive and 
negative twist as well as monotonic twist profiles, as have been 
more commonly studied. We have shown that it is not possible 
to determine the onset of instability by any simple criterion such 
as achieving a critical twist or critical helicity. Nevertheless, Fig. 
|4] suggests that the location of the instability is coincident with 
the peak twist (the largest absolute twist). 

A distribution of heating events has been calculated, both 
for the case of a single loop which undergoes repeated stressing 
and relaxation, and for an ensemble of identical loops which are 
randomly stressed (and for intermediate cases). For a sufficiently 
large number of events, a statistically-stable distribution of event 
sizes is obtained. As expected, the smallest events are the most 
common, and there is (as noted by Browning & Van der Linden 
2003) a minimum event size for given loop parameters, perhaps 
corresponding to an "elemental nanoflare". More surprisingly, 
there is a second peak of event frequency at intermediate magni- 
tudes, although this is somewhat reduced in size if an ensemble 
of loops is considered. This can be explained as follows. The 
first peak (at minimum energy) occurs simply because the range 
of energies near the minimum naturally encompasses the largest 
part of the instability threshold curve (because of the flatness of 
the minimum). The second peak is found because the instabil- 
ity threshold is most likely to be crossed in the part near to the 
region of constant-o;. 
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An ensemble of identical loops - in which each individual 
loop starts from a randomly chosen initial state and is stressed 
until it undergoes a single relaxation event - is investigated in 
more detail. The distribution of heating events is qualitatively 
similar, although, the secondary peak of event frequency is much 
lower, and the decay of frequency with increasing energy is 
much flatter. A significant result is that, although the distribution 
is not a simple power law, the high-energy part of the distribution 
is well approximated by a power law with an index of around - 
8.3, which is considerably steeper than the minimum required 
for nanoflare heating to be eff'ective (-2). 

The model requires that the field is sometimes unstable - al- 
though most of the time, the field profile will be well within the 
stable region. Typically, the dimensional a value is of magni- 
tude 1-2, leading to dimensional values of o- * 1 - 2 Mm ' (for 
a loop radius of 1 Mm). Note that this is the maximum value 
which could be found, and usually we would expect lower val- 
ues. This is consistent with observations; for example Regnier & 
Priest (2007) find a magnitudes around 1 Mm"'. Furthermore, 
a consequence of the fact that the fields are predicted to fluc- 
tuate between stable and unstable states is that we predict a 
value for the average horizontal field component (on average, 
this will be somewhat less than the value at marginal stability). 
At threshold, we find BJB~ is around 0.5, which means that its 
average value should be around 0.25. This agrees very well with 
the limit on this quantity required by Parker (1988), in order for 
the Poynting flux from the photosphere to match coronal heat- 
ing requirements. Hence, we have an explanation for the critical 
shear value predicted by Parker This also implies that the av- 
erage heating flux from our model will be sufficient for coronal 
heating. Indeed, the fluxes derived from our results agree (espe- 
cially if we raise the correlation time for photospheric motions) 
with the required values. 

Relaxation theory makes no prediction about the spatial dis- 
tribution of energy dissipation. However, this can be determined 
from numerical simulations. Recently, Hood et al. (2009) have 
shown that heating is well-distributed across the loop volume, as 
the current sheet, associated with the nonlinear kink instability, 
stretches and fragments, thereby filling the loop cross section. 
This has implications for the observed emission. 

A further important consideration is that we study a single 
loop of fixed dimensions, with repeated stressing and relaxation. 
In practice, coronal heating events will occur in regions of field 
of varying sizes. A large flare will inevitably involve a large mag- 
netic field volume, whereas nanoflares may involve small sets of 
fieldlines. This could be accounted for in our model by allowing 
the loop radius to be randomly distributed also; the effect would 
be to convolve a set of our energy distributions. A second limita- 
tion of the calculation is that we have assumed a uniform random 
walk in a space. This is not realistic, since photospheric twisting 
motions are likely to be correlated across the loop cross section 
and therefore changes in which the twist in the outer later is sim- 
ilar to that in the core are much more likely: in other words, there 
should be a positive correlation between changes in a\ (the core 
current) and changes in ao (the current in the outer layer), rather 
than these being independent random variables as assumed here. 
This will be considered further in future work. 

As well as the improvements mentioned above, the model 
will be enhanced such that there is, more realistically, zero net- 
current carried by the loop. At present, in most cases the po- 
tential layer outside the loop contains azimuthal field. A current 
neutralisation layer will be inserted between the loop and the en- 
velope - external magnetic fields will have an axial direction only 



(Hood et al. 2009), representing a response to localised photo- 
spheric twist motions. 

The model has included a number of simplifications, and 
at this stage is more of a "proof of principle" showing that 
a distribution of heating events can be produced from an 
(almost) ab initio coronal heating model. These initial results 
demonstrate the viability of the model for further research. 
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Appendix A: Expressions for loop properties 

Expressions for some key quantities ((^>, K and W) are given here. 
For compactness, these are given only for o-i # and ai i^ 0, while spe- 
cial cases (e.g., ai =0) must be dealt with separately. Expressions for 
constant-a fields can be recovered by setting o-i =0-2^ which gives more 
familiar formulae. The superscripts and subscripts that accompany each 
quantity term denote the upper and lower radial bounds over which the 
quantity is calculated. 



A.1. 


Average magnetic twist 
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RiJi{\a,\R{) 






(^)2 



0-2L[Fo(|Q'2l^l)-/^o(|Q'2l^2)J 

\R2Fd\a2\R2) - RiFi{\a2\Ri}] 



2cr2LR2[B2J,(\a2\R2) + C2yi(|f^2l^2)]log(|) 
B2Fo{\a2\R2)\Ri - R^] 



(A.1) 



(A.2) 



(A.3) 



A.2. Magnetic helicity 



K^' = o-i^^\RlJ^(\ai\Ri) + R'^r^(\a,\Ri) 



lo'il 



„ IcTnnLB} 
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A.3. Magnetic energy 
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Appendix B: Magnetic field profiles for a selection 
of cy-space points 
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Fig. B.l. The magnetic field profiles, B- (solid) and Eg (dashed), for 
the six Q--space points identified in the top left panel of Fig.|4l 



